DEMO: HF calculations on small organic molecules by PySCF#
We need to install a few things before we get started
%pip install -q pyscf py3DMol plotly kaleido
Note: you may need to restart the kernel to use updated packages.
from pyscf import gto # Used to define a molecule
from pyscf import scf # Used to perform HF calculations
from pyscf import mp # Used to perform Møller–Plesset PT calculations
from pyscf import cc # Used to perfrom Coupled Cluster calculations
from pyscf import mcscf # Used to perform multireference calculations
Build and Visualize Molecules#
mol_xyz = """H 1.2194 -0.1652 2.1600
C 0.6825 -0.0924 1.2087
C -0.7075 -0.0352 1.1973
H -1.2644 -0.0630 2.1393
C -1.3898 0.0572 -0.0114
H -2.4836 0.1021 -0.0204
C -0.6824 0.0925 -1.2088
H -1.2194 0.1652 -2.1599
C 0.7075 0.0352 -1.1973
H 1.2641 0.0628 -2.1395
C 1.3899 -0.0572 0.0114
H 2.4836 -0.1022 0.0205"""
mol = gto.M(atom=mol_xyz, basis="sto3g", verbose=3)
import py3Dmol
xyz_view = py3Dmol.view(width=400,height=400)
xyz_view.addModel(mol.tostring(format="xyz"),'xyz')
xyz_view.setStyle({'stick':{}, "sphere":{"radius":0.4}})
xyz_view.setBackgroundColor('0xeeeeee')
xyz_view.show()
#
# Use your mouse to interact with the molecule!
#
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
Run calculations#
# Run HF and save the results
mf = scf.RHF(mol).run()
# In Jupyter notebooks you can hover over methods to see docstrings or you can print the docstrings out!
#print(mf.__doc__)
converged SCF energy = -227.890481261003
Visualize Results#
import plotly.express as px
# Plot the MO Occupations
fig = px.line(y=mf.mo_occ, markers=True, title="Molecular Orbital (MO) Occupations")
fig.update_layout(xaxis_title="Orbital Index (0-Based)", yaxis_title="MO Occupation")
fig.show()
# Plot the MO Energies (i.e. eigenvalues of the Fock matrix)
fig = px.line(y=mf.mo_energy, markers=True, title="Molecular Orbital (MO) Energies")
fig.update_layout(xaxis_title="Orbital Index (0-Based)", yaxis_title="MO Energies")
fig.show()
from pyscf.tools import cubegen
cubegen.orbital(mol, 'mo.cube', mf.mo_coeff[:,18]); # 42 electrons so 21 occupied orbitals
cube_view = py3Dmol.view(width=400,height=400)
cube_view.addVolumetricData(open("mo.cube").read(), "cube", {'isoval': -0.03, 'color': "red", 'opacity': 0.85})
cube_view.addVolumetricData(open("mo.cube").read(), "cube", {'isoval': 0.03, 'color': "blue", 'opacity': 0.85})
cube_view.addModel(mol.tostring(format="xyz"),'xyz')
cube_view.setStyle({'stick':{}, "sphere":{"radius":0.4}})
cube_view.setBackgroundColor('0xeeeeee')
cube_view.show()
#
# Use your mouse to interact with the molecule!
#
3Dmol.js failed to load for some reason. Please check your browser console for error messages.